#Author: Ilana M Ventura					
#Date: September 2024															  
#Step 7: Make figures from SAS output

library(foreign)
library(tidyverse)
library(mclust)
library(dplyr)
library(readstata13)
library(longCatEDA)
library(magick)
library(webshot)
library(kableExtra)
library(knitr)
library(haven)
library(githubinstall)
library(ggpattern)
library(ggpol)
library(ggpubr)
library(ggrepel)
library(ggsignif)


##Figure S1

plot.data <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/table1.sas7bdat')

levels(as.factor(plot.data$race))
plot.data$race<-(factor(plot.data$race, levels = c("White*", "Hisp","Black", "Asian",
                                                   "MENA*", "SOR*",  "AmIndAN","NHoPI" , "Multiracial")))

plot.data <-subset(plot.data, race !="Multiracial")

plot.data$Treatment<-(factor(plot.data$Treatment, levels = c("Separate", "Combined")))


ggplot(plot.data, aes(fill=Treatment, y=Est, x=race)) + 
  geom_bar(position="dodge", stat="identity")+
  geom_errorbar(aes(ymin=Est, ymax=Est+stderr), width=.1, position=position_dodge(.9)) +
  theme_bw() +
  #scale_fill_grey(start = .3, end = .8)+
  xlab(' ')+ 
  ylim(0, .8)+
  ylab("Proportion")+
  ggsave("H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/figs/figure_S1.png", 
         width = 6, height = 4)



#fig 2a
  
  plot.data2a <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig2_a.sas7bdat')
  
  levels(as.factor(plot.data2a$Group))
  plot.data2a$Group<-(factor(plot.data2a$Group, levels = c("AmIndAN", "SOR","Black", "Hisp", "White"   )))
  

  ggplot(data=plot.data2a ,   aes(x =Var2 , y = Mean, ymin = LowerCLMean, ymax = UpperCLMean ))+
    geom_pointrange(aes(col= Group), position = position_dodge(.7), size=.2)+
    theme_bw() +
    coord_flip()+
    facet_wrap(~Var1, scales = "free") +
    theme(legend.position="bottom", 
          legend.title=element_blank(),
          strip.placement = "outside",
          #strip.text.y = element_text(face = "bold"), 
          panel.grid.major.y = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA), 
          strip.background = element_blank())+
    geom_hline(yintercept=0, linetype='longdash')+
    xlab('')+ 
    ylab("Mean Difference, 95% CI")+
    scale_color_hue(direction = -1, h.start=0)+
    guides(color = guide_legend(reverse=TRUE))+
  ggsave("H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/figs/figure2_a.png", 
         width = 3.5, height = 2)
  
  
  
  #fig 2b
  
  plot.data2b <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig2_b.sas7bdat')
  
  
  levels(as.factor(plot.data2b$Var2))
  plot.data2b$Var2<-(factor(plot.data2b$Var2, levels = c(  "Less than $30K" , "$30K-$59K" , "$60K-$99K" ,
                                                         "$100K or more" , "English",  "Bilingual" , "Spanish" , 
                                                         "Democrat" , "Independent",   "Republican" ,  "US Born" ,  "Born in LatAm/Sp",
                                                         "Less than HS", "High School",  "College", 
                                                         "Light",  "Medium" ,"Dark", "18-29" , "30-44" , "45-64" ,  "65+" )))
  
  levels(as.factor(plot.data2b$Var1))
  plot.data2b$Var1<-(factor(plot.data2b$Var1, levels = c( "Skin Color Scale" , "Age Group" ,  "Education",
                                                          "Reported Income" ,  "Partisanship" ,    "Language" ,  "Nativity")))
  
  levels(as.factor(plot.data2b$Group))
  plot.data2b$Group<-(factor(plot.data2b$Group, levels = c("AmIndAN", "SOR","Black", "Hisp", "White"   )))
  
  ggplot(data=plot.data2b ,   aes(x =Var2 , y = Mean, ymin = LowerCLMean, ymax = UpperCLMean ))+
    geom_pointrange(aes(col= Group), position = position_dodge(.7), size=.2)+
    theme_bw() +
    coord_flip()+
    facet_wrap(~Var1, scales = "free") +
    # scale_colour_grey(start = 0.8, end = 0)+
    theme(legend.position="bottom", 
          legend.title=element_blank(),
          strip.placement = "outside",
          #strip.text.y = element_text(face = "bold"), 
          panel.grid.major.y = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA), 
          strip.background = element_blank())+
    geom_hline(yintercept=0, linetype='longdash')+
    xlab('')+ 
    ylab("Mean Difference, 95% CI")+
    scale_color_hue(direction = -1, h.start=0)+
    guides(color = guide_legend(reverse=TRUE))+
  ggsave("H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/figs/figure2_b.png", 
         width = 7, height = 6)
  
  

  
  
#fig 3 hisp
  
plot.data3h <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig3_hisp.sas7bdat')
  
plot.data3h$Group<-(factor(plot.data3h$Group, levels = c("Separate", "Combined")))
  
levels(as.factor(plot.data3h$Variable))
plot.data3h$Variable<-(factor(plot.data3h$Variable, levels = c( "White",  "Black" , "HispAlone", "SOR", "MENA" , "NHoPI" ,"AmIndAN" , "Asian"  )))
  
  
hisp<- ggplot(plot.data3h, aes(fill=Group, y=Mean, x=Variable)) + 
    geom_bar(position="dodge", stat="identity", width=.7)+
    geom_errorbar(aes(ymin=Mean, ymax=Mean+StdErr), width=.1, position=position_dodge(.9)) +
    theme_bw() +
    xlab(' ')+ 
    ylim(0, 1)+
    ylab("")+
   geom_signif(y_position = c(.68,.1,.72,.34,.05,.05,.15,.05),
               xmin = c(.8,1.8,2.8, 3.8,4.8,5.8,6.8,7.8), 
               xmax = c(1.2,2.2,3.2,4.2,5.2,6.2,7.2,8.2),  
              annotations=c("*","+","*","*","*","ns","*","+"),  tip_length = .01, vjust=-.1)+
  scale_x_discrete(guide = guide_axis(angle = 30)) +
  guides(fill=guide_legend(title="Treatment"))+
  ggtitle("Hispanic")
hisp
  

#fig 5 black

plot.data3b <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig3_black.sas7bdat')

plot.data3b$Group<-(factor(plot.data3b$Group, levels = c("Separate", "Combined")))

levels(as.factor(plot.data3b$Variable))
plot.data3b$Variable<-(factor(plot.data3b$Variable, levels = c( "White",  "BlackAlone" , "Hisp", "SOR", "MENA" , "NHoPI" ,"AmIndAN" , "Asian"  )))


black<- ggplot(plot.data3b, aes(fill=Group, y=Mean, x=Variable)) + 
  geom_bar(position="dodge", stat="identity", width=.7)+
  geom_errorbar(aes(ymin=Mean, ymax=Mean+StdErr), width=.1, position=position_dodge(.9)) +
  theme_bw() +
  xlab(' ')+ 
  ylim(0, 1)+
  ylab("")+
  geom_signif(y_position = c(.1,.95,.1,.05,.05,.05,.08,.05),
              xmin = c(.8,1.8,2.8, 3.8,4.8,5.8,6.8,7.8), 
              xmax = c(1.2,2.2,3.2,4.2,5.2,6.2,7.2,8.2),  
              annotations=c("ns","ns","ns","ns","ns","ns","ns","ns"),  tip_length = .01, vjust=-.1)+
  scale_x_discrete(guide = guide_axis(angle = 30))+
  guides(fill=guide_legend(title="Treatment"))+
  ggtitle("Black")
black




#fig 5 White

plot.data3w <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig3_white.sas7bdat')

plot.data3w$Group<-(factor(plot.data3w$Group, levels = c("Separate", "Combined")))

levels(as.factor(plot.data3w$Variable))
plot.data3w$Variable<-(factor(plot.data3w$Variable, levels = c( "WhiteAlone",  "Black" , "Hisp", "SOR", "MENA" , "NHoPI" ,"AmIndAN" , "Asian"  )))


White<- ggplot(plot.data3w, aes(fill=Group, y=Mean, x=Variable)) + 
  geom_bar(position="dodge", stat="identity", width=.7)+
  geom_errorbar(aes(ymin=Mean, ymax=Mean+StdErr), width=.1, position=position_dodge(.9)) +
  theme_bw() +
  xlab(' ')+ 
  ylim(0, 1)+
  ylab("")+
  geom_signif(y_position = c(.92,.05,.2,.05,.05,.05,.08,.05),
              xmin = c(.8,1.8,2.8, 3.8,4.8,5.8,6.8,7.8), 
              xmax = c(1.2,2.2,3.2,4.2,5.2,6.2,7.2,8.2),  
              annotations=c("*","ns","*","ns","*","ns","*","ns"),  tip_length = .01, vjust=-.1)+
  scale_x_discrete(guide = guide_axis(angle = 30))+
  guides(fill=guide_legend(title="Treatment"))+
  ggtitle("White")
White


#fig 5 SOR

plot.data3s <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig3_SOR.sas7bdat')

plot.data3s$Group<-(factor(plot.data3s$Group, levels = c("Separate", "Combined")))

levels(as.factor(plot.data3s$Variable))
plot.data3s$Variable<-(factor(plot.data3s$Variable, levels = c( "White",  "Black" , "Hisp", "SORAlone", "MENA" , "NHoPI" ,"AmIndAN" , "Asian"  )))


SOR<- ggplot(plot.data3s, aes(fill=Group, y=Mean, x=Variable)) + 
  geom_bar(position="dodge", stat="identity", width=.7)+
  geom_errorbar(aes(ymin=Mean, ymax=Mean+StdErr), width=.1, position=position_dodge(.9)) +
  theme_bw() +
  xlab(' ')+ 
  ylim(0, 1)+
  ylab("")+
  geom_signif(y_position = c(.3,.15,.95,.8,.15,.05,.1,.05),
              xmin = c(.8,1.8,2.8, 3.8,4.8,5.8,6.8,7.8), 
              xmax = c(1.2,2.2,3.2,4.2,5.2,6.2,7.2,8.2),  
              annotations=c("*","*","*","*","*","ns","ns","ns"),  tip_length = .01, vjust=-.1)+
  scale_x_discrete(guide = guide_axis(angle = 30))+
  guides(fill=guide_legend(title="Treatment"))+
  ggtitle("SOR")
SOR


#fig 5 AmInAN

plot.data3i <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/fig3_amindan.sas7bdat')

plot.data3i$Group<-(factor(plot.data3i$Group, levels = c("Separate", "Combined")))

levels(as.factor(plot.data3i$Variable))
plot.data3i$Variable<-(factor(plot.data3i$Variable, levels = c( "White",  "Black" , "Hisp", "SOR", "MENA" , "NHoPI" ,"AmIndANAlone" , "Asian"  )))


AmInAN<- ggplot(plot.data3i, aes(fill=Group, y=Mean, x=Variable)) + 
  geom_bar(position="dodge", stat="identity", width=.7)+
  geom_errorbar(aes(ymin=Mean, ymax=Mean+StdErr), width=.1, position=position_dodge(.9)) +
  theme_bw() +
  xlab(' ')+ 
  ylim(0, 1)+
  ylab("")+
  geom_signif(y_position = c(.7,.2,.5,.1,.05,.05,.3,.1),
              xmin = c(.8,1.8,2.8, 3.8,4.8,5.8,6.8,7.8), 
              xmax = c(1.2,2.2,3.2,4.2,5.2,6.2,7.2,8.2),  
              annotations=c("*","ns","*","ns","ns","ns","ns","ns"),  tip_length = .01, vjust=-.1)+
  scale_x_discrete(guide = guide_axis(angle = 30))+
  guides(fill=guide_legend(title="Treatment"))+
  ggtitle("AmInAN")
AmInAN
  

  
  figure3b <- ggarrange(hisp, White, SOR, black,AmInAN,
                       ncol = 2, nrow = 3,  common.legend = TRUE)
  
  figure3b+
    ggsave("H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/figs/figure3.png", 
           width = 8, height = 9)

  
######interactions pred probs (Fig S2 and S3)
  
pred_int <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/pred_interactions.sas7bdat')

levels(as.factor(pred_int$type))
pred_int$type<-(factor(pred_int$type, levels = c("Separate", "Combined")))

levels(as.factor(pred_int$pred))
pred_int$pred<-(factor(pred_int$pred, levels = c( "Origin" ,"Nativity"  , "Discrimination",  "Skin color"  )))

pred_int <- pred_int %>% filter(pred!="Discrimination")

ggplot(data=pred_int ,   aes(x =int , y = est_response, ymin = L, ymax = U ))+
  geom_pointrange(aes(col= type), position = position_dodge(.7), size=.2)+
  theme_bw() +
  coord_flip()+
  facet_wrap(~pred, scales = "free") +
  theme(legend.position="bottom", 
        legend.title=element_blank(),
        strip.placement = "outside",
        #strip.text.y = element_text(face = "bold"), 
        panel.grid.major.y = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA), 
        strip.background = element_blank(), 
        axis.title.x = element_text(vjust=-2))+
  xlab('')+ 
  ylim(0, .3)+
  ylab("Predicted Probability, 95% CI")+
  ggsave("H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/figs/figure_S2.png", 
         width = 6, height = 3)

###Indig hisp


pred_int_ind <- read_sas('H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/pred_interactionsind.sas7bdat')

levels(as.factor(pred_int_ind$type))
pred_int_ind$type<-(factor(pred_int_ind$type, levels = c("Separate", "Combined")))

levels(as.factor(pred_int_ind$pred))
pred_int_ind$pred<-(factor(pred_int_ind$pred, levels = c( "Origin" ,"Nativity"  , "Discrimination",  "Skin color"  )))

pred_int_ind <- pred_int_ind %>% filter(pred!="Discrimination")

ggplot(data=pred_int_ind ,   aes(x =int , y = est_response, ymin = L, ymax = U ))+
  geom_pointrange(aes(col= type), position = position_dodge(.7), size=.2)+
  theme_bw() +
  coord_flip()+
  facet_wrap(~pred, scales = "free") +
  theme(legend.position="bottom", 
        legend.title=element_blank(),
        strip.placement = "outside",
        #strip.text.y = element_text(face = "bold"), 
        panel.grid.major.y = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA), 
        strip.background = element_blank(), 
        axis.title.x = element_text(vjust=-2))+
  xlab('')+ 
  ylim(0, .3)+
  ylab("Predicted Probability, 95% CI")+
  ggsave("H:/UChicago/Census Experiment/Census Experiment Replication/Replication Files/out/figs/figure_S3.png", 
         width = 6, height = 3)

